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Abstract. We describe a new image co-addition tool, AWAIC, to support 
the creation of a digital Image Atlas from the multiple frame exposures acquired 
with the Wide-field Infrared Survey Explorer (WISE). AWAIC includes prepara- 
tory steps such as frame background matching and outlier detection using robust 
frame-stack statistics. Frame co-addition is based on using the detector's Point 
Response Function (PRF) as an interpolation kernel. This kernel reduces the im- 
pact of prior-masked pixels; enables the creation of an optimal matched filtered 
product for point source detection; and most important, it allows for resolution 
enhancement (HiRes) to yield a model of the sky that is consistent with the 
observations to within measurement error. The HiRes functionality allows for 
non-isoplanatic PRFs, prior noise-variance weighting, uncertainty estimation, 
and includes a ringing-suppression algorithm. AWAIC also supports the popu- 
lar overlap-area weighted interpolation method, and is generic enough for use 
on any astronomical image data that supports the FITS and WCS standards. 



1. Introduction 

The goal of image co-addition is to optimally combine a set of (usually dithered) 
exposures to create an accurate representation of the sky, given that all instru- 
mental signatures, glitches, and cosmic-rays have been properly removed. By 
"optimally" , we mean a method which maximizes the signal-to-noise ratio (SNR) 
given prior knowledge of the statistical distribution of the input measurements. 

The Wide-field Infrared Survey Explorer (WISE) mission will be generating 
over a million exposures (or frames) over the sky. WISE is a NASA Midex 
mission scheduled for launch in late 2009. It will survey the entire sky at 3.3, 
4.7, 12, and 23;um with sensitivities up to three orders of magnitude beyond 
those achieved with previous all-sky surveys. For details on the scientific goals, 
requirements, instrument and mission design, see Mainzer et al. 2005. One of 
the primary products from WISE is a digital Image Atlas that combines the 
multiple 8.8 second, 47' x 47' frame exposures within predefined tiles over the 
sky. To support this, we have developed a suite of software modules collectively 
referred to as AWAIC for execution in the automated processing pipeline at the 
Infrared Processing and Analysis Center. The modules are written in ANSI- 
compliant C and wrapped into a Perl script, and will be made portable in the 
near future. 

Here we review AWAIC's co-addition algorithms, products, and extension to 
resolution enhancement (HiRes). It's important to note that HiRes is not in the 
WISE baseline plan. It was implemented primarily to support offline research. 
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INPUTS : (images in FITS format) 

- instrumentally/astrometrically calibrated image frames 

- bad-pixel masks 

- uncertainty (sigma) images from prior noise model 



Throughput matching (multiplicative) to scale input 
frames to a common photometric zero-point 



Frame background-level matching (additive) 

1 

Outlier detection & masking 

Interpolate frames onto a common grid 
using overlap-area weighting method 



Final Product Generator. Deliver 
to public through IRSA at IPAC 



OUTPUTS : 

- main co-add intensity images 

- depth-of-coverage maps 

- uncertainty images 

- all 4096 X 4096; si .3757pixel 



Quality Assurance: backgrounds, noise, 
coverage/outlier stats, stats for 
uncertainty validation 



Threshold in each interpolated pixel stack 
using robust statistics. Masks updated 



Co-addition of all unmasked pixels using AWAIC: 

- fast re-projection and distortion correction 

- interpolation: PRF weighted averaging with 
optional inverse variance weighting 



Figure 1. Processing steps in WISE frame co-addition pipeline. 



The statistical robustness and performance of algorithms will be addressed in 
more detail in future papers. 



2. WISE Frame Co-addition Pipeline 

Figure [1] gives an overview of the co- addition steps. It is assumed that the input 
science frames have been preprocessed to remove instrumental signatures and 
their pointing refined in some WCS using an astrometric catalog. Accompany- 
ing bad-pixel masks (in 32-bit integer format) and prior-uncertainty frames are 
optional. The frames are assumed to overlap with some predefined footprint (or 
tile) on the sky. This also defines the dimensions of the co-add products. The 
uncertainty frames store 1-cr values for each pixel. These are expected to be 
initiated upstream from a noise model specific to the detector and then propa- 
gated and updated as the instrumental calibrations are applied. The uncertain- 
ties are used for optional inverse- variance weighting of the input measurements, 
and for computing co-add flux uncertainties. If bad-pixel masks are specified, 
a bit-string template is used to select which conditions to flag against. The 
corresponding pixels in the science frames are then omitted from co-addition. 

The first (optional) step is to scale the frame pixel values to a common 
photometric zero-point using calibration zero-point information in each FITS 
header. Currently, the software reads a zero-point in magnitudes stored in the 
"MAGZP" keyword. The common (or target) zero point is then written to the 
FITS headers of the co-add products to enable the calibration of photometric 
measurements. Frame overlap (or background-level) matching and outlier de- 
tection is then performed. These are described in § [3]] Since these initial steps 
modify the input frame and mask pixel values, local copies of the frames and 
masks are made to avoid overwriting the originals. After outliers have been 
flagged in the input masks, the frames are ready for co-addition. All "good" 
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(unmasked) pixels are reprojected and interpolated onto an upsampled output 
grid. Details are described in § [J]] The reprojection uses a fast input-to-output 
plane coordinate transformation that implicitly corrects for focal plane distor- 
tion if represented in the input FITS headers. The Simple Imaging Polynomial 
(SIP) convention for distortion is assumed (Shupe et al. 2005). 

The outputs from AWAIC are the main intensity image, a depth-of-coverage 
map, a 1-a uncertainty image based on input priors, an image of the outlier 
locations, and optionally if the overlap-area interpolation method was used, an 
image of the data-derived uncertainty computed from the standard deviation in 
each interpolated pixel stack and appropriately scaled by the depth-of-coverage. 
AWAIC also produces a wealth of Quality Assurance (QA) metrics and plots 
over pre-specified regions of the co-add footprint. These include background 
noise estimates, coverage and outlier statistics, and metrics to validate co-add 
flux uncertainties using tests. 

3. Preparatory Steps 

3.1. Background Matching 

Frame exposures taken at different times usually show variations in background 
levels due to, for example, instrumental transients, changing environments, scat- 
tered and stray light. The goal is to obtain seamless (and/or smooth) transitions 
between frames near their overlap regions prior to co-addition. We will want to 
equalize background levels from frame to frame, but at the same time preserve 
natural background variations and structures as much as possible. We have 
implemented the following simple method: 

1. Fit a robust plane to each frame. By "robust", we mean immune to the 
presence of bright sources and extended structure. Our goal is to cap- 
ture the global underlying background level in a frame. There are of 
course cases where structure may span over most of a frame, and hence 
the background will be over-represented. The planar fit is parameterized 
by 2 = f{x, y), where z is the background level, and x, y are frame-pixel 
coordinates. 

2. The robust planar fits are subtracted from each respective frame. This 
effectively flattens the frames and places them on a zero-baseline. 

3. Compute either: (i) the global median M of all frame pixels contributing 
to the co-add footprint, or, (ii) a median plane from all the planar fits. 
The latter attempts to capture any natural background gradient over the 
co-add footprint. 

4. Add this global median M (a constant) to each of the "zero-level" frames 
(from 2), or, in the case of a median plane, extend it over the co-add 
region and then add it self-consistently to each input frame. The frames 
have now been matched to a common background level (or gradient). This 
will be more or less representative of the natural background over the co- 
add footprint region. 

This method ensures continuity of the background across the footprint re- 
gion after co-addition. This not only improves a co-add's esthetic appearance (by 
minimizing frame level offsets between overlaps) , but also makes it self-consistent 
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for scientific purposes. It's important to note that instrumental transients or 
improper instrumental calibration (e.g., bad flat-fielding) can also manifest as 
gradients across frames. Therefore, one needs to be sure that any retained global 
gradient is purely astrophysical. 

The above also includes a method to ameliorate biases from the possible 
presence of bright extended structure. Presence of extended structure (e.g., a 
galaxy) over a frame is first searched for by thresholding on the ratio of quantile 
differences in the pixel distribution, e.g., Qd = [go.84 — Qo.5]/[qo.5 — ^o.ie]- Values 
of Qd 2 usually indicate a highly skewed distribution and hence contamination 
from bright extended structure. If detected, a frame is partitioned into a grid 
and only those regions with the lowest median background value are used to 
perform the robust planar fitting. This method still has its limitations, but it 
extends the robustness of the algorithm. 

3.2. Outlier Detection and Masking 

The goal of outlier detection is to identify frame pixel measurements of the same 
location on the sky which appear inconsistent with the (bulk) remainder of the 
sample at that location. This assumes multiple frame exposures of the same 
region of sky are available. Potential outliers include cosmic rays, latents (im- 
age persistence) , instrumental artifacts (including bad pixels) , poorly registered 
frames from gross pointing errors, supernovae, asteroids, and basically anything 
that has moved or varied appreciably with respect to the inertial sky over the 
observation span of a set of overlapping frames. 

In summary, the method involves first projecting and interpolating each 
input frame onto a common grid with user-specified pixel scale optimized for 
the detector's Point Spread Function (PSF) size. The interpolation is performed 
using the overlap-area weighting method (analogous to using a top hat kernel). 
This accentuates and localizes the outliers for optimal detection (e.g., cosmic 
ray spikes). When all frames have been interpolated, robust estimates of the 
first and second moments are computed for each interpolated pixel stack j. We 
adopt the sample median (med) , and the Median Absolute Deviation (MAD) as 
a proxy for the dispersion: 



where pi is the value of the r interpolated pixel within stack j. The factor of 
1.4826 is the correction necessary for consistency with the standard deviation of 
a Normal distribution in the large sample limit. The MAD estimator is relatively 
immune to the presence of outliers where it exhibits a breakdown point of 50%, 
i.e., more than half the measurements in a sample will need to be declared 
outliers before the MAD gives an arbitrarily large error. 

The final step involves re-projecting and rc-intcrpolating each input pixel 
again, but now testing each for outlier status against other values in its stack 
using the pre-computed robust metrics. A pixel with value pi is declared an 
outlier if for given "upper" (uthres) and "lower" {Ithres) t^^il thresholds, either of 
the following is satisfied: 



(Tj = 1.4826 med — med{pi}\} 




Pi > med{pi] + UthresCTj 
Pi < med{pi] - khresCTj 



(2) 
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Figure 2. One-dimensional schematic of stacking method for detecting out- 
hers for well sampled (left) and under-sampled (right) cases. The input pixel 
marked " X " contains signal from a source and is in danger of being flagged 
when outlier detection is performed at location j in the output grid. 



If declared an outlier, a bit is set in the accompanying frame mask for use 
downstream. The algorithm also includes an adaptive thresholding method in 
that if a pixel is likely to contain "real" signal (e.g., from a source), the upper 
threshold is automatically inflated by a specified amount to avoid (or reduce the 
incidence of) outlier flagging at that location. To distinguish between what's 
real or not, we generate a background subtracted medianSNH co-add using all 
the input pixels. The background and local noise are computed using spatial 
median filtering and quantile differencing: a ~ go.5 — Qo.w respectively. The idea 
here is that since these metrics are relatively outlier resistant, a large median 
pixel value in the co-add (or SNR derived therefrom) is likely to contain signal 
associated with a source. Therefore, when flagging outliers using Eq. [21 we also 
threshold on the SNR co-add to determine if Uthres should be inflated. 

We require typically at least five samples (overlapping pixels) in a stack for 
the above method to be reliable. This is because the MAD measure for a, even 
though robust, can itself be noisy when the sample size is small. Simulations 
show that for the MAD to achieve the same accuracy as the most optimal esti- 
mator of a for a normal distribution (i.e., the sample standard deviation), the 
sample size needs be ~ 2.74x larger. A noisy a will adversely affect the ability 
to perform reliable outlier detection. Another requirement to ensure good reli- 
ability is to have good sampling of the instrumental PSF, i.e., at the Nyquist 
rate or better. When well sampled, more detector pixels in a stack can be made 
to align within the span of the PSF, and any pixel variations due to PSF shape 
are minimized. On the other hand, a PSF which is grossly under-sampled can 
artificially increase the scatter in a stack, with the consequence of erroneously 
flagging pixels containing true source signal. Figure [2] illustrates these concepts. 
The WISE detectors are all slightly better than critically sampled. Simulations 
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have shown that for depths-of-coverage of eight or more, (where eight is the 
median depth for WISE when scanning across the echptic), we expect to detect 
outhers to completeness and rehabihty levels of ^ 80% for a nominal threshold 

of 5(7. 



4. Co-addition using PRF Interpolation 

One of the interpolation methods in AWAIC involves using the detector's Point 
Response Function (PRF) as the interpolation kernel. The PRF is simply the 
instrumental PSF convolved with the pixel response. When knowledge of the 
intra-pixel responsivity is absent, the pixel response is assumed to be uniform, 
i.e., a top hat. The PRF is what one usually measures off an image using the 
profiles of point sources. Each pixel can be thought as collecting light from its 
vicinity with an efficiency described by the PRF. 

The PRF can be used to estimate the flux at any point in space as follows. 
In general, the flux in an output pixel j is estimated by combining the input 
detector pixel measurements Di using PRF-weighted averaging: 

where rij is the value of the PRF from input pixel i at the location of output 
pixel J. The r^j are volume normalized to unity, i.e., for each = 1. This 

will ensure flux is conserved. The inverse- variance weights (Xjof) are optional 
and default to 1. The Oi can be fed into AWAIC as 1-cj uncertainty frames, e.g., 
as propagated from a prior noise model. The sums in Eq. [3] are over all input 
pixels in all input frames. With multiple overlapping input frames, this will 
result in a co-add. The l-cr uncertainty in the co-add pixel flux /j, as derived 
from Eq. [3] is given by 

1/2 



E2 2 



(4) 



where Wij = {'rij/o'f)/Yl,i "^ijl^l- Equation[l]assumes the measurement errors (in 
the Di) are uncorrelated. Note that this represents the co-add flux uncertainty 
based on priors. With N j overlapping input frames and assuming Oi = constant 
throughout, it's not difficult to show that Eq. U] scales as: aj ~ ai/ y/NjPj, 
where Pj = 1/J2i''"ij is a characteristic of the detector's PRF, usually referred 
to as the effective number of "noise pixels". This scaling also assumes that the 
PRF is isoplanatic (has fixed shape over the focal plane) so that Pj = constant. 
Furthermore, the depth-of-coverage at co-add pixel j is given by the sum of all 
overlapping PRF contributions at that location: Nj = This effectively 

indicates how many times a point on the sky was visited by the PRF of a "good" 
detector pixel i, i.e., not rejected due to prior-masking. If no input pixels were 
masked, this reduces to the number of frame overlaps, Nf. 

In general, the PRF is usually non-isoplanatic, especially for large detector 
arrays. AWAIC allows for a list of spatially varying PRFs to be specified, where 
each PRF corresponds to some pre-determined region (e.g., a partition of a 
square grid) on the detector focal plane. 
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Figure 3. Schematic of PRF interpolation for a single input pixel. 

Equation [3] can be compared to the popular pixel overlap-area weighting 
method, e.g., as implemented in the Montagqjtool. In fact if the PSF is grossly 
under-sampled, then the PRF is effectively a top hat spanning one detector pixel. 
The interpolation as described above then reduces to overlap-area weighted av- 
eraging where the interpolation weights rij become the input(i)-to-output(j) 
pixel overlap areas aij. Incidentally, AWAIC also implements the overlap-area 
weighting method, in case detector PRFs are not available. 

Figure [3] shows a schematic of a detector PRF mapped onto the co-add 
output grid. The PRF boundary is shown as the dashed circle and is centered on 
the detector pixel. To ensure accurate mapping of PRF pixels and interpolation 
onto the co-add grid, a finer cell-grid composed of "pixel cells" is set up internally. 
The cell size can be selected according to the accuracy to which the PRF can be 
positioned. The PRF is subject to thermal fluctuations in the optical system as 
well as pointing errors if multiple frames are being combined. Therefore, it does 
not make sense to have a cell- grid finer than the measured positional accuracy of 
the PRF. The PRF pixels are mapped onto the cell-grid frame by first projecting 
the center of the detector pixel with distortion correction if necessary, and then 
using a local transformation with rotation to determine the positions of the PRF 
pixels in the cell-grid. The value of a PRF- weighted detector pixel flux rijDi in 
a co-add cell pixel j is then computed using either a nearest-neighbor match, 
or, the overlap-area weighting method. The latter is more accurate but slower. 
After all the input pixels with their PRFs have been mapped, the internal co-add 
cells are down-sampled to the desired output co-add pixel sizes. 

There are three advantages to using the PRF as an interpolation kernel. 
First, it reduces the impact of masked (missing) input pixels if the data are well 
sampled, even close to Nyquist. This is because the PRF tails of neighboring 
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"good" pixels can overlap and stretch into the bad pixel locations to effectively 
give a non-zero coverage and signal there in the co-add. Second, Eqs [3] and H] can 
be used to define a linear matched filter optimized for point source detection. 
This effectively represents a cross-covariance of a point source template (the 
PRF) with the input data. It leads to smoothing of the high-frequency noise 
without affecting the point source signal sought. In other words, the SNR per 
pixel in the co-add is maximized for detecting point source peaks. The inclusion 
of inverse-variance weighting further ensures that the SNR is maximized since 
it implies the co-add fluxes will be maximum likelihood estimates for normally 
distributed data. The creation of co-adds which are also optimized for source- 
detection will benefit projects (e.g., WISE) where a source catalog is also a 
release product. The third advantage is that the PRF allows for resolution 
enhancement by "deconvolving" its effects from the input data. 

Use of the PRF as an interpolation kernel also has its pitfalls, at least for the 
process of co-add generation. The operation defined by Eq. [3] leads to a "smooth- 
ing" of the input data in the co-add grid. This smoothing is minimized for a 
top-hat PRF spanning one detector pixel (equivalent to overlap-area weighting) . 
This leads to smearing of the input pixel signals and one consequence is that 
cosmic rays can masquerade as point sources (albeit with narrower width) if not 
properly masked. For point sources with Gaussian profiles, their effective width 
will increase by a factor of ~ \/2. Furthermore, a broad kernel will cause the 
noise to be spatially correlated in the co-add, typically on scales (correlation 
lengths) approaching the PRF size. Correlations are minimized for top-hat ker- 
nels. Both the effects of flux smearing and correlated noise must be accounted 
for in photometric measurements off the co-add, both in profile fitting and aper- 
ture photometry. The compensation for flux smearing can be handled through 
an appropriate aperture correction. Ignorance of correlated noise will cause pho- 
tometric uncertainties to be underestimated. Methods on how to account for 
correlated noise in photometry will be discussed in a future paper. 

5. Extension to Resolution Enhancement 

We now describe a generic framework for co-addition with optional resolution 
enhancement (HiRes). Above we referred to the concept of combining frames 
to create a co-add. The HiRes problem asks the reverse: what model or rep- 
resentation of the sky propagates through the measurement process to yield 
the observations within measurement error? As a reminder, the measurement 
process is a filtering operation performed by the instrument's PRF: 

sky (truth) (g) PSF (g □ sampling measurements. (5) 

PRF 

Our goal is to infer a plausible model of the sky or "truth" given the instrumental 
effects. 

5.1. The Maximum Correlation Method 

The HiRes algorithm in AWAIC is based on the Maximum Correlation Method 
(MCM). This was originally implemented to boost the scientific return of data 
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from IRAS approximately 20 years ago (Aumann et al. 1990; Fowler & Au- 
mann 1994), and is still provided as an online service to users. We have now 
implemented MCM in a form which is suitable for use on any imaging data that 
are compatible with the FITS and WCS standards, and the SIP convention for 
distortion. The versatility of MCM is that it implicitly generates, as its very 
first step (or first iteration), a PRF-interpolated co-add as described above. The 
algorithm is as follows. 

1. First we begin with a flat model image of ones, i.e., a "maximally corre- 
lated" image: 

/;=o = ivj, (6) 

where the subscript j refers to a pixel in the upsampled output grid, and 
n refers to the iteration number. This starting image is a first guess at the 
"truth" that we plan to reconstruct. Obviously this is a bad approxima- 
tion, since it represents what we know without any measurements having 
been used yet. We could instead have used prior information as the start- 
ing model if it was available. 

2. Next, we use the detector PRF(s) to "observe" this model image, or predict 
the input detector measurements. Starting with n = 1, the predicted flux 
in each detector pixel i is obtained by a "convolution": 

J 

where is the response (PRF value) of pixel i at the location of output 
model pixel j. Eq. [7]is a tensor inner product of the model image with the 
flipped PRF (see below for why we need to flip the PRF). It may not be 
a true convolution since the kernel rij may be non-isoplanatic. 

3. Correction factors are computed for each detector pixel i by dividing their 
measured flux, Di, by those predicted from the model (Eq. [7|): 

= (8) 

4. For each model pixel j, all "contributing" correction factors, i.e., con- 
tributed by the overlapping PRFs of all neighboring detector pixels i are 
averaged using response- weighted averaging (with optional I /erf weight- 
ing): 

5. Finally, the model image pixels are multiplied by their respective averaged 
correction factors (Eq. [9]) to obtain new refined estimates of the model 
fluxes: 

/7 = /r'^7- (10) 

If we are after a simple PRF-interpolated co-add, we terminate the process 
at step 5. In fact, Eq. [9] is analogous to the co-addition equation (Eq. [3|) in 
that a starting model image with /j^ = 1 implies a correction factor Kl = Di 

since a PRF volume-normalized to unity predicts = 1 (Eq. [7|). Therefore 
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after the first (n = 1) iteration of MCM, co-add fluxes will automatically result: 



If we desire resolution enhancement, the above process is iterated, where the 
model image from step 5 is used to re-predict the measurements in step 2. This 
process of iteratively refining the model continues until the model reproduces the 
measurements to within the noise, i.e., the predictions from Eq. [7]are consistent 
with the measurements Di. If input prior uncertainties (di) are available, this 
convergence can be formally checked using a global t^st that uses all the 
input detector pixels: 



where we expect Xn — ^ (the number of degrees of freedom, = the number 
of input pixels). Alternatively, convergence can also be checked by examining 
the correction factors for each detector pixel (Eq. [8|) , where we expect K^^ ~ 1 
within the noise, or, via the averaged correction factors (Eq. llOp . where C" 1 
after many iterations. An image of the latter can be generated by the software 
at each iteration. Iterating much further beyond the initial signs of convergence 
has the potential of introducing unnecessary (and usually unaesthetic) detail in 
the model. This is important to ensure a parsimonious HiRes solution. 

Therefore, it is an algorithmic property of MCM that it only modifies (or 
de-correlates) a flat starting model image to the extent necessary to make it 
reproduce the measurements within the noise. A PRF-interpolated co-add (from 
the first MCM iteration) will generally not satisfy the measurements after it is 
"convolved" with the detector PRFs, i.e., when subject to the measurement 
process (Eq. [5]). 

As a detail, the input PRFs are first flipped in x and y (or equivalently 
rotated by 180°) when HiRes'ing is performed (n > 1). This is to conform 
to the usual rules of convolution and assumes the input PRFs were made by 
combining images of point sources observed with the same detector in the same 
native x-y pixel frame. For PRF-interpolated co-adds however (that terminate 
at n = 1), the PRFs are not flipped since a cross-covariance with the input 
data is instead needed. The PRFs here are used as matched filters to generate 
products optimized for point source detection (see §113). 

It is also worth noting that MCM reduces to the classic Richardson-Lucy 
(RL) method if the following are assumed: (i) the PRE is isoplanatic so that a 
constant kernel allows for Fourier-based deconvolution methods to be used; (ii) 
the inverse- variance weighting of measurement correction factors is disabled from 
the PRE- weighted averaging (Eq. ED, or equivalently if all the input variances 
erf are assumed equal. This implies the solution will converge to the maximum 
likelihood estimate for data that are Poisson distributed. With inverse-variance 
weighting included, the solution converges to the maximum likelihood estimate 
for Gaussian distributed data. This is usually always satisfied for astronomical 
image data in the limit of high photon counts; (iii) there is no explicit testing for 
global convergence at each iteration by checking, for example, that the solution 
reproduces the data within measurement error (Eqs [7] and [TT|) . This criterion 
was indeed suggested by Lucy (1974), although it is seldom used in modern 
implementations of the RL method. 
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In the absence of prior information for the starting model, MCM impUcitly 
gives a solution which is the "smoothest" possible, i.e., has maximal entropy. 
This should be compared to maximum entropy methods (e.g., Cornwell & Evans 
1985) which attempt to minimize the of the differences between the data and 
the convolved model, with an additional constraint imposing smoothness of the 
solution. MCM requires no explicit smoothness constraint. MCM can indeed use 
a regularizing constraint in the form of non-flat starting model, (e.g., an image 
of the sky from another detector or wavelength), but this jettisons the idea of an 
image with maximally correlated pixels, and the refined model image will not 
be the smoothest possible. Smoothness is important because it can be used to 
convey fidelity in a model. In general, the solution to a deconvolution problem is 
not unique, especially in the presence of noise. Many models can be made to fit 
the data, and many methods invoke regularization techniques in order to select 
a plausible solution. A consequence is that some methods give more structure 
or detail than necessary to satisfy the data, and there is no guarantee that this 
structure is genuine. MCM adopts the Occam's razor approach. Given no prior 
constraints (apart from satisfying the input data), MCM will always converge 
on the simplest solution. This will be the smoothest possible. 

5.2. The CFV Diagnostic 

A powerful diagnostic from MCM is the Correction Factor Variance (CFV). 
This represents the variance about the PRF-weighted average correction factor 
(Eq. [9]) at a location in the output grid for iteration n: VJ^ = {Kf)j — {Ki)"], or 



where Wij = {r-ij / af) / ^.^rij / af , and the detector-pixel correction factors Kf 
were defined in Eq. [51 At early iterations, the CFV is generally high everywhere 
because spatial structure has not yet been resolved, and the model contradicts 
the measurements when subject to the measurement process. If after conver- 
gence, all the detector-pixel measurements contributing a non-zero response at 
some location j agreed exactly with their predicted fluxes (Eq. [7]) , then all the 
Kf would be ~ 1 and the CFV (V^") at that location would be zero. Areas 
with a relatively large CFV indicate the presence of input pixel measurements 
which do not agree with the majority of the other measurements (e.g., outliers). 
It could also indicate noisy data, saturated data, regions where the FRF is not 
a good match (e.g., erroneously broad), or that a fleld has not yet converged 
and would benefit from further iteration. By thresholding the CFV, one can 
therefore create a mask for a HiRes image to assist in photometry, e.g., to avoid 
outliers and unreliable detections from amplified noise fluctuations (see below). 

Example CFVs for the M51 galaxy are shown in Figure [6j The correspond- 
ing co-add and HiRes'd images appear in Figure O To illustrate the above 
concepts, outlying input measurements were not masked in the left and middle 
images of Figure [H These refer to iteration levels n = 1 and n = 40 respec- 
tively, with the latter corresponding to convergence. The CFV image on the far 
right was created from data with outliers detected and masked a priori using 
the algorithm described in § 13.2.1 
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Apart from providing a qualitative diagnostic, the CFV can also be used 
to compute a posteriori (data-derived) uncertainties in the pixel fluxes /" in a 
HiRes image. In general, the 1-a uncertainty at iteration n can be written in 
terms of the CFV as: 

= ^"/ry^TVE^' (13) 

where the sum is over the responses from all measurements i at output pixel 
j, i.e., the effective depth-of-coverage. c" is a correction factor to account for 
re-distribution of noise power across spatial frequencies from one iteration to 
the next. At low iterations, power is relatively high at low frequencies, i.e., the 
noise is correlated across pixels. As iterations increase, noise is de-correlated 
and power migrates to high frequencies. The spectrum approaches that of white 
noise, provided the input measurement noise was spectrally white. For n = 1 
(giving a co-add), = Xj ^fPj^ where Pj is the effective number of noise pixels 
defined in §2]] With written this way, Eq.ll3lbecomes equivalent to the co-add 
pixel uncertainty defined in Eq. HI In general, the at any iteration n > 1 can 
be approximated from the output image products as: 

C^iJAfs[/"] ^^^^ 



where (Jrms is the standard-deviation (or some robust equivalent) of the pixel 
noise fluctuations within a "source-free" stationary background region with 
~uniform depth-of-coverage in the /" image. The denominator is the mean 
(or median) of Eq. [13] with c" = 1 in the same region. At the time of writing, 
AWAIC only computes an image of (t"[c" = 1], since it can be quite subjective 
on how the source-free stationary background is chosen. If such doesn't exist, 
background fitting may be required with urms computed from the fit residuals. 
The user can then rescale the a^[c^ = 1] image using the estimate of c" from 
Eq. [TH This will give pixel uncertainties which are more or less statistically 
compatible with noise fluctuations in the HiRes'd image. Pixel SNRs will also 
be the maximum possible since MCM would have converged to the maximum 
likelihood estimate for data that were Gaussian distributed. With the correct 
value of c", a user then has an estimate of the flux uncertainty anywhere in the 
HiRes'd image, including at the location of sources. This will allow one to esti- 
mate uncertainties in source photometry. Noise correlations are also expected to 
be minimal in a converged HiRes image, or negligible if products were created 
with ringing suppression turned on (see below). 

5.3. Ringing Suppression 

Like most deconvolution methods, MCM can lead to ringing artifacts in the 
model image. This limits super-resolution, i.e., when attempting to go well 
beyond the diffraction limit of an imaging system. In general, ringing occurs 
because the reconstruction process tries to make the model image agree with 
the "true" scene with access to only the low spatial frequency components com- 
prising the data. The input data are usually band- limited, and information 
beyond some high spatial frequency cutoff can never be recovered. The best we 
can ever reconstruct is a "low-pass filtered" version of the truth, with the filter 
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Input image frames 

Each frame : compute background using blocl< median 
or mode filter over N x N grid w/ Gaussian smoothing 



Mosaic the background 
frames 



- subtract background frames from input frames 

- reset all negative pixel values to zero 

- add small positive offset 6-1 O ^" 

{=> positivity constraint with a "flux bias") 



Using a flat model (starting) image, HIRes 
hckgnd-SUbtraGted frames until convergence 



Add background mosaic 
to HiRes products 



Using new HiRes as starting image and 
oriainal frames (wl bckands). HiRes aaain 
to 4-10 iterations to restore (de-bias) noise 
structure and reproduce measurements 








Output HiRes'd mosaic with ringing 
suppressed. Residual ringing depends 
on complexity of background 






Figure 4. Ringing suppression algorithm. 



determined by the maximum spatial frequency the observations provide. This 
includes the finite samphng by pixels. A hard high frequency cutoff will lead to 
sinc-like oscillations in real image space. The magnitude of the ringing depends 
on the strength of a source relative to the local background intensity level. 

It is no accident that a solution with ringing is the smoothest (and sim- 
plest) solution possible with MCM. Anything smoother (with more low frequency 
power) will not satisfy the measurements when subject to the measurement pro- 
cess (Eq. [S]). However, since a large number of less-smooth solutions can re- 
produce the observations, those without ringing are generally more desirable. 
Therefore, we relax our request for the smoothest image and use prior knowl- 
edge that the background and (desired) source fluxes are physically distinct and 
separable. There have been numerous approaches that have used this philosophy 
(e.g., Lucy 1994). In brief, the ringing suppression algorithm in AWAIC first 
generates an image of the slowly varying background for each input frame on 
some specific scale using median filtering; this is subtracted from the respective 
input frames to create the "source" images; negative noise fluctuations are set 
to zero, and a tiny positive offset added; MCM is then run on the background- 
subtracted images until convergence; the background images are combined and 
then added to the HiRes'd source-image product. This operation enforces a 
positivity constraint for reconstruction of the source signals. It ensures that 
source flux won't ring against an essentially zero background level so power can 
be forced into high spatial frequencies. After the background has been added 
to the HiRes'd source-only product, MCM is re-executed for several iterations 
using this as the starting model image and the original frames as input. This 
step re-adjusts the solution and attempts to restore the intrinsic noise properties 
of the HiRes process, i.e., what one would have obtained if no background were 
removed or positivity constraint enforced. It ensures that photometric uncer- 
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Figure 5. M51 from Spitzer-MIPS 24/Lim. Left: co-add (after 1 iteration); 
Middle: HiRes after 10 iterations; Right: HiRes after 40 iterations. 

tainties don't become biased and the final solution adequately reproduces the 
measurements within the noise. Figure [4] gives an overview of the above steps. 

5.4. HiRes in Practice 

Like most deconvolution methods, MCM does not alter the information content 
of the input image data. The signal and noise at a given frequency are scaled 
approximately together, keeping the SNR fixed. The process just re-emphasizes 
different parts of the frequency spectrum to make images more amenable to a 
certain kind of examination, e.g., for detecting previously unresolved objects 
and thereby increasing the completeness of surveys. 

For optimal HiRes'ing, the input data will have to adequately sample the 
instrumental PSF to at least better than the Nyquist sampling frequency 2i'c, 
where Vc is the maximum frequency cutoff inherent in the PSF. For a simple 
diffraction-limited system with aperture diameter D, Uc (x D/X and corresponds 
to the full width at half maximum (FWHM) of an Airy beam. Even if the 
detector pixels undersample the PSF (below Nyquist), redundant coverage with 
N randomly dithered frames can help recover the high spatial frequencies, since 
the average sampling will scale as ~ 1/A^ of an input pixel. The better the 
sampling, the better the HiRes algorithm is at improving spatial resolution. 
For imaging data from the Spitzer IRAC and MIPS detectors with typically 
SNR^ 5/pixel and 10 frame overlaps, our HiRes algorithm reduces the FWHM 
of the effective PRF to ~ 0.35X/D - a factor of almost 3 below the diffraction 
limit. This corresponds to almost an order of magnitude increase in flux per solid 
angle for a Gaussian profile. This enhancement assumes accurate knowledge of 
the PRF over the focal plane. 

An example output from AWAIC at three MCM iteration levels is shown in 
Figure O At high iterations, point source ringing starts to appear. The ringing 
around the satellite dwarf galaxy at the bottom is aggravated because the core 
is saturated in the data, and the PRF used for HiRes'ing (which is derived from 
unsaturated data) is not a good match. "Flat" core profiles in the data, due 
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Figure 6. Correction Factor Variance (CFV) images for M51 whose inten- 
sity images were shown in Fig. [5j For a description of the CFV, see § 15.2.1 
Left: CFV after 1 iteration with outlying measurements purposefully retained; 
Middle: CFV after 40 iterations with outlying measurements also purposefully 
retained; Right: CFV after 40 iterations but with outliers masked (omitted) 
prior to HiRes'ing. Darkest regions correspond to lowest values of the CFV 
{Vj < 0.1), and the brightest to highest values {Vj > 100). 

to either saturation or improperly corrected non-linearity, will contain relatively 
more power in the side-lobes than the actual PRF. When this PRF is used for 
HiRes'ing, these side-lobes will manifest as ringing artifacts in the HiRes image 
in order for it to reproduce the observations on "convolution" with the PRF. 
Even though the ringing suppression algorithm was turned on in this example, 
ringing is still seen around other point sources. This is because these sources 
are superimposed on the extended structure of the galaxy. This structure acts 
like an elevated background against which point sources can ring. The ringing 
suppression algorithm relies on accurate estimation of the local background, and 
this can be difficult when complex structure is involved, as it is here. 



6. Summary and Future Work 

We have given a broad overview of the algorithms implemented in a new generic 
co-addition/HiRes'ing tool. The goal is to produce high fidelity science qual- 
ity products with uncertainty estimates and metrics for validation thereof. The 
HiRes (MCM) algorithm contains considerable improvement over previous meth- 
ods in that it includes a posteriori uncertainty estimation, statistically motivated 
convergence criteria, a powerful diagnostic (the CFV) to locate inconsistencies 
in the input data and assess the overall quality of HiRes solutions, and the abil- 
ity to handle non-isoplanatic PRFs. Algorithms will be discussed in more detail 
in future papers. Future work will explore methods to accelerate convergence 
in MCM, the ability to handle time-dependent PRFs (e.g., adapted to variable 
seeing), and an analysis of the completeness, reliability, and photometric accu- 
racy of sources detected in HiRes'd images, especially in confused fields. More 
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examples, analyses, user-interface details, and animations of MCM can be found 
at http : //web . ipac . caltech. edu/staff/f mas ci/home/wise/awaic .html 
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